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ABSTRACT 



Aims. To study the Reynolds stresses which describe turbulent momentum transport from turbulence affected by large-scale shear 
and rotation. 

Methods. Three-dimensional numerical simulations are used to study turbulent transport under the influences of large-scale shear and 
rotation in homogeneous, isotropically forced turbulence. We study three cases: one with only shear, and two others where in addition 
to shear, rotation is present. These cases differ by the angle (0 or 90°) the rotation vector makes with respect to the z-direction. Two 
subsets of runs are performed with both values of 6 where either rotation or shear is kept constant. When only shear is present, the 
off-diagonal stress can be described by turbulent viscosity whereas if the system also rotates, nondiffusive contributions (A-effect) to 
the stress can arise. Comparison of the direct simulations are made with analytical results from a simple closure model. 
Results. We find that the turbulent viscosity is of the order of the first order smoothing result in the parameter regime studied and that 
for sufficiently large Reynolds numbers the Strouhal number, describing the ratio of correlation to turnover times, is roughly 1.5. This 
is consistent with the closure model based on the minimal tau-approximation which produces a reasonable fit to the simulation data 
for similar Strouhal numbers. In the cases where rotation is present, separating the diffusive and nondiffusive components of the stress 
turns out to be challenging but taking the results at face value, we can obtain nondiffusive contributions of the order of 0.1 times the 
turbulent viscosity. We also find that the simple closure model is able to reproduce most of the qualitative features of the numerical 
results provided that the Strouhal number is of the order of unity. 
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1. Introduction 

■ Turbulent angular momentum transport is considered to be of 
importance in various astrophysical objects, such as accretion 
disks (e.g. Balbus & Hawley 1998 ) and convectively unstable 
layers within stars (e.g. Riidiger 1989) where they take part in 

] shaping the internal rotation profile of the object. Due to the 
immense numerical requirements, direct global simulations of 
these systems are not yet feasible in large quantities, although 
some simulations are able to capture many features of, for exam- 
ple, the solar rotation profile (see e.g. Robinson & Chan 120011 

• Brun & ToomreEOOll Miesch et al. l2006l lIOOSl l. However, in 
many cases it would be desirable to be able to use simplified, and 
computationally less demanding, mean-field models where the 
effects of small-scale turbulence are accurately parameterized 
in some collective way. Although a rich literature of mean-field 
models of solar internal rotation exist (e.g. Brandenburg et al. 
[19921 Kiiker et al. 1993 Kitchatinov & Rudiger|2005l Rempel 
120051 ). many of the models use simple and often untested param- 
eterizations of the turbulent quantities. 

Parameterizing turbulence entails a closure model for the tur- 
bulent correlations. One of the most often used closures in astro- 
physics is the a-prescription (Shakura & Sunyaev 1973 ), widely 
used in accretion disk theory, which relates the turbulent viscos- 
ity to the local gas pressure. This very simple parameterization 
allows analytically tractable solutions of accretion disk structure 
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but suffers from the drawback that important physics, such as 
magnetic fields, are not taken into account. More recently, dy- 
namical turbulent closure models yielding all the relevant com- 
ponents of the Reynolds and Maxwell stresses in simplified se- 
tups with homogeneous magnetized shear flows have started to 
appear (e.g. Ogilvie i2003i Pessah et al. 120061 ). 

For these more sophisticated models to be useful, they need 
to be validated somehow. The obvious validation method is 
to compare to numerical simulations that work in a parameter 
regime which is as similar as possible. Comparisons of closure 
models with numerical simulations of the magnetorotational in- 
stability in periodic local slab geometry have appeared recently 
(e.g. Pessah et al. ,2008. Liljestrom et al. 2009) and yield en- 
couraging results in that the closure models are able to capture, 
at least qualitatively, many of the main features of the numer- 
ical simulations. However, these studies concentrate on a sit- 
uation where the turbulence is generated via the magnetorota- 
tional inst ability (Velikhov 1959; Chandrasekhar 119611 Balbus 
& Hawley 119911 ). the nonlinear behaviour of which is still not 
very well understood (e.g. Fromang et al. 120071 ). 

In the present paper we avoid these complications and con- 
sider only the simplest possible hydrodynamic case and assume 
that isotropic and homogeneous background turbulence akeady 
exists in the system upon which large-scale shear and rotation 
can be imposed. Such turbulence can be generated by using 
suitable forcing in the Navier-Stokes equation. Although flows 
of this kind are not likely to occur in nature, they are perfect 
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testbeds for the developement and testing of turbulent closure 
models. 

The present paper is a continuation to an earlier study 
(Kapyla & Brandenburg 2008 hereafter KB08; see also Kapyla 
& Brandenburg .2007 J where simulations of anisotropic homo- 
geneous turbulence under the influence of rotation were com- 
pared to a simple closure model applying the so-called tau- 
approximation (hereafter MTA, see e.g. Blackman & Field 
125021 120031 Brandenburg et al. l2004l l. In the present study we 
adopt an isotropic forcing and add a large-scale shear flow using 
the shearing box approximation. Such a setup allows us to deter- 
mine the turbulent viscosity (see preliminary results in Kapyla & 
Brandenburg 2002 and Kapyla et al. 2009). The imposed shear 
flow also introduces anisotropy into the turbulence which, under 
the influence of rotation, can produce additional non-diffusive 
Reynolds stresses (Leprovost & Kim 2007 2008a 20081^, 
which are more commonly known as the A-effect (Krause & 
Rudiger 1974. Rudiger [T980|1989l ). We make an effort to sep- 
arate the diffusive and nondiffusive contributions from the nu- 
merical data. As in KB08, one of the main goals of the study 
is to compare the simulation results to a simple closure model, 
similar to that introduced by Ogilvie (2003 ). 

The remainder of the paper is organised as follows: the mod- 
els used in the study are presented in Sect. [H whereas Sects. [3] 
and|4|give the results and conclusions of the study, respectively. 



2. The model and methods 

2. 1. Basic equations 

We model compressible hydrodynamic turbulence in shearing 
periodic cube of size {^tt)^. The gas obeys an isothermal equa- 
tion of state characterized by a constant speed of sound, Cg . The 
continuity and Navier-Stokes equations read 



V\np 
Vt 



= -v-u. 



(1) 



vu 

— = -SU^y - c2 V In p - 2 X [/ + /vise + /force, (2) 

where V/Vt = d/dt + {U + TJq) ■ V denotes the advective 
derivative and Uq — (0, Sx, 0) is the imposed shear flow. U is 
the velocity, p is the density, /vise is the viscous force, and /foreo 
is the forcing function. Compressibility is retained but low Mach 
number flows, i.e. Ma = Urms/cs ~ 0.05 — 0.1, where Wrms 
is the average root mean square velocity, are considered. The 
viscous force is given by 



/vise = i^(v2[/+iVV-t/ + 2S- Vlnp 



where i/ is the kinematic viscosity and 



(3) 



(4) 



is the traceless rate of strain tensor The forcing function /force 
is given by 



fix, t) - Re{iV/fc(t) exp[ifc(t) • x - i(j){t)]} 



(5) 



where x is the position vector, = foCsikcs/StY^"^ is a nor- 
malization factor, /o is the forcing amplitude, k = \k\, 6t is the 



length of the time step, and — tt < < tt a random delta- 
correlated phase. The vector /^ is given by 



fk 



k X e 



y/k^ - [k ■ e)2 



(6) 



where e is an arbitrary unit vector. Thus, /;- describes nonhelical 
transversal waves with |/fcp = 1, where k is chosen randomly 
from a predefined range in the vicinity of the average wavenum- 
ber fcf/fci = 5 at each time step. Here fci is the wavenumber 
corresponding to the domain size, and fcf is the wavenumber of 
the energy-carrying scale. The choice fcf/fci = 5 is somewhat 
arbitrary but provides a clear scale separation between the tur- 
bulent eddies and the system size. 

The numerical simulations were performed with the PENCIL 
which uses sixth-order accurate finite differences 
in space, and a third-order accurate time-stepping scheme 
(Brandenburg & Dobler 2002; Brandenburg |2003"l l . Resolutions 
up to 1024^ grid points were used. 

2.2. Dimensionless units and parameters 
We obtain nondimensional variables by setting 

Cs = fci = po = 1- (7) 
This means that the units of length, time, and density are 

[x]=k-\ [t]^{c,k,)-\ [p]=po. (8) 

However, in what follows we present the results in explicitly di- 
mensionless form using the quantities above. 

The strengths of shear, rotation, and viscous effects are mea- 
sured by the shear, Coriolis and Reynolds numbers, respectively, 
based on the forcing scale as 



Sh 



S 



^rmsfcf 



2 17q ^rms 

Co — , Re 



^rms^f 



l^kf 



(9) 



See Fig.[T]for a typical snapshot from a high resolution run with 
Re « 387. 

2.3. Coordinate system, averaging, and error estimates 

The simulated domain can be thought to represent a small rect- 
angular portion of a spherical body of gas. We choose (x, y, z) to 
correspond to (0, 0, r) of spherical coordinates. With this choice 
the rotation vector can be written as 



rj = ilo(— sin0, 0, cos ( 



(10) 



where Q is the angle between the rotation axis and the local ver- 
tical (z-) direction, i.e. the colatitude. We consider two cases: 
Q — Q and Q — 90°. We can consider these cases to rep- 
resent the north pole and the equator of a rotating star, re- 
spectively. However, if an additional a shear flow of the form 
TJq = (0, Sx, 0) is introduced into the system there are multi- 
ple ways to interpret the physical system that is described by the 
model. 

Firstly, the case Q — Q can be considered to describe a local 
portion of a disk rotating around a central object sufficiently far 
away that the effects of curvature can be neglected. Then the 
rotation profile of the disk is characterized by f2 cx i?^^ where R 
is the radius and q = —S/flo- Now, a Keplerian rotation profile 
is obtained for g = 1.5, a flat profile, such as those observed in 
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t= 50.0 




Fig. 1. Velocity component Uy, in the units of the sound speed, 
from the Run B6 with Re = 387, Sh = -0.20, and Co = 0, 
resolution 1024^. 

many galaxies, is given by q = 1, and for perfectly rigid rotation 
we have g = 0. 

Alternatively, the shear flow can be understood to represent 
either radial or latitudinal shear in a convection zone of a star 
This approach has been used by Leprovost & Kim (I2007l l who 

consider that in the case 6 = 0, C/^"'' corresponds to latitudinal 
shear near the equator, and in the case 9 = 90° to the radial shear 
in near the pole of the star 

Since the turbulence is homogeneous, volume averages are 
employed and denoted by overbars. An additional time aver- 
age over the statistically saturated state of the simulation is also 
taken. Errors are estimated by dividing the time series into three 
equally long parts and computing mean values for each part in- 
dividually. The largest departure from the mean value computed 
for the whole time series is taken to represent the error 

2.4. Reynolds stresses from the minimal 
tau-approximation 

We follow here the same procedure as in KB08 and derive a 
simple analytical model for the Reynolds stresses in the case 
of homogeneous turbulence under the influences of rotation and 
shear. Although we present the model in terms of the mini- 
mal tau-approximation (see, e.g. Blackman & Field |200"2l 12003 1 
Brandenburg et al. 2004; KB08) the closure used here is quite 
similar to that originally presented by Ogilvie (120031 see also 
Garaud & Ogilvie |20051 l. A more detailed comparison is given 
below. 

One of the main purposes of this study is to compare the 
results from the closure model with numerical simulations. Since 
the numerical setup is homogeneous, it is sufficient to compare 
the volume averaged data with a closure model with no spatial 
extent. In this case, the Navier-Stokes equations yield 

il, = ~SUJ,y ~ U.jdjU, ~ eld, Inp - 2eakniUk + (11) 

where the dot represents time derivative, and /; describes both 
the viscous force and external forcing. Now we decompose the 
velocity as Ui = Ui + Ui, where Ui is the average velocity and 



Ui the fluctuation. It is straightforward to derive the equation for 
iii from Eq. (fTTT i which allows us to derive an evolution equation 
for the Reynolds stress Qij = UjUj, 

Qij — SSyiQxj — SSyjQxi '^^ilk^lQki ^ "^^jlk^lQkj 

-c^Uidj \np + Ujdi \np+QijV ■ u + fiUj + fjUi, (12) 

where the overbars denote averaging. In the minimal tau- 
approximation closure scheme the nonlinear terms in Eq. ST% 
are modeled collectively by a relaxation term 

— c^Uidj \np + Ujdi Inp + Qij'V ■ u = —r^^Qij, (13) 

where r is a relaxation time. Terms involving the forcing and 
viscosity, fi, can be parameterized in the same spirit with 

7^ + l^^ = rf'Q["\ (14) 

where q'"' is the equilibrium solution in the absence of shear 
and rotation. Throughout this paper we assume that the time 
scale associated with the forcing is equal to the relaxation time, 
i.e. T = Tf . We note that there is no compelling theoretical argu- 
ment to enforce this equality but we rather use it for the sake of 
simplicity. Thus we arrive at the equation 

Qij SSyiQxj SSyjQxi '^^ilk^lQki '^^jlk^lQkj 

-T-\Q,,-Qf/). (15) 

We use the same values of and S as in the simulations and take 
the quantity q'"' from a nonrotating run with Co = Sh = 0. The 
free parameter in the model is the relaxation time, r, which can 
be represented in terms of the Strouhal number 

St = TUrnish, (16) 

which is the ratio of correlation and turnover times. 

2.4.1 . Comparison to the Ogilvie (2003) model 

In the model of Ogilvie ( I2003I I the Reynolds stresses in the hy- 
drodynamical case are given by 

Qij + linear terms = 

CiVgL-ig,, + C2VQL-\Q^3 - ^QS^J), (i7) 

where d are dimensionless parameters of the order of unity, Q 
is the trace of the Reynolds tensor and i is a lenght scale, e.g. the 
system size. The first term on the rhs can be identified as a relax- 
ation term similar to that used in the minimal tau-approximation 

where we have used ^/Q = Urms, and L = 27r/fci. Comparing 
with Eq. (O, we see that St = 27rfcf (Cifci)-i = lOTrCf ^ for 
our value of k[/ki = 5. 

Our Eq. (fTSI l also lacks the term C2VQL^^{Qij ^ ^QSij) 
which appears in the model of Ogilvie. This term describes 
isotropization of the turbulence, and a similar term can be added 
to the minimal tau-approximation in the form Tj^^ (Qy — QSij), 
where tj could have a value unequal to t. However, after exper- 
imenting with such a term, we did not find substantially better 
agreement with simulation results. Thus, in order to keep our 
model as simple as possible with the minimum amount of free 
parameters, also this term is neglected in our analysis. 
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Table 1. Summary of the different sets of simulations. The val- 
ues of Re, Co, and Sh are given in terms of Mi,„s from a run with 
Co = Sh = 0. 



Set 


Re 


Co 


-Sh 


e 


q 




A 5. 


, . . 345 





0.1 








B 8. 


, . . 387 





^ 0.2 








C 


=s 27 





0.01 . . .0.24 








D 


^ 75 





0.01... 0.23 








E 


=s 24 


-1.27. . . 1.28 


fa 0.16 





-5... 


1.9 


F 


=s 24 


f» 0.32 


-0.29... 0.64 





-5... 


1.9 


G 


=s 24 


-1.27. . . 1.28 


^ 0.16 


90° 


-5... 


1.9 


H 


=s 24 


f» 0.32 


-0.29... 0.64 


90° 


-5... 


1.9 



Eq. ( |20l i. This formulation is supported by numerical results al- 
though analytical studies often yield somewhat different results 
(e.g. Yousef et al. 2003 and references therein). 

The correlation time can be related to the turnover time of 
the turbulence via the Strouhal number, Eq. (fT6b . If we assume 
that St 1, as is suggested by numerical turbulence models 
(e.g. Brandenburg & Subramanian 120051 |2007| l similar to ours, 
the turbulent viscosity is given by 

I^tO = Iwrmsfcf"^ (21) 

In what follows we use lyto as the normalization for the turbulent 
viscosity. However, later on we allow St 7^ 1 and write vt — 
Stfto in order to obtain an independent estimate of St from the 
numerical simulations. 



We note that if the saturated solutions of Eqs. ( fTSl l and ( fTTI i 
are independent of time, the two models yield the same results 
provided that C2 is zero. However, if Q 7^ in the saturated 
state, e.g. if the solution is oscillatory, then differences will occur 
because the relaxation time in our model is based on a constant 
Uims and not on the Q that emerges as a results of time integra- 
tion of the model itself. We also stress that in the original work 
of Ogilvie (I2003I I, the turbulence due to the shear flow, rotation, 
and magnetic fields was studied. In our case, however, the turbu- 
lence pre-exists due to the external forcing and we mainly study 
the effects of shear and rotation on this background turbulence. 

3. Results 

We study the Reynolds stresses from three different systems: one 
where only shear flow is present and two others where nonzero 
rotation is present with either 9 = Qoid ^ 90°. For the first sys- 
tem we perform four different sets of simulations (see Table[T]for 
a summary of the different sets of runs) where we either fix Sh 
(Sets A and B) or Re (Sets C and D). In the setups with rotation 
we perform two sets of calculations for both values of 6 where 
we fix either Sh (Sets E and G) or Co (Sets F and H). In Sets E 
to H we vary the parameter q in the range —5 < q < 1.9. More 
detailed descriptions of the runs can be found in Tables IbTI to 
IB. 31 Furthermore, we compute the stresses from corresponding 
closure models making use of the stationary solution of Eq. (flST l 
given in Appendix lAl 

3. 1. Case il = 0: turbulent viscosity 

Consider first the case where rotation is absent and the only 
large-scale flow is the imposed shear U — (0, xS, 0) with 
isotropic background turbulence. The Reynolds stress generated 
by the shear can be represented by the expression 

Qij = -vt{Ui^j + U j^i) = -vtSdixSjy, (19) 

where i^t is the turbulent viscosity. The expression iT% can be 
considered to be valid for weak shear On the other hand, we can 
estimate the turbulent viscosity from 

z/to = (20) 

where t is the correlation time and is the average turbulent 
velocity squared. Equation ( |20l ) is the same as the one that can be 
derived for the magnetic diffusivity using first order smoothing 
approximation (e.g. Krause & Radler |19801 l. Here we have as- 
sumed a turbulent magnetic Prandtl number of unity to arrive at 



3.1 .1 . Simulation results 

We find that in the non-rotating case it is difficult to avoid large- 
scale vorticity generation (see also Kapyla et al. 2009). This phe- 
nomenon is likely to be related to a vorticity dynamo discussed 
in the analytical studies of Elperin et al. ( 120031 12007I I. Similar 
large-scale structures have been observed earlier in the numeri- 
cal works of Yousef et al. (2008a' '2008b') using an independent 
method. The problem becomes increasingly worse as the shear 
is increased. Thus we need to limit the range of Sh and carefully 
analyze the data so that the effect of the vorticity dynamo on the 
turbulence is minimized (see the related discussion in Mitra et 
al. |2009l l. What this means is that we limit the time averaging to 
a range where turbulence is established but where the large-scale 
flow is still weak in comparison to the turbulent rms-velocity. We 
stress that the vorticity dynamo occurs only if turbulence already 
exists in the system, i.e. it does not arise if the forcing is turned 
off in the simulations. 

The only nonzero off-diagonal component in the simulations 
is now Qxy that can be interpreted in terms of turbulent viscos- 
ity. This is also consistent with symmetry arguments. Once a 
suitable averaging interval has been found, the turbulent viscos- 
ity is obtained from Eq. (T% . If the first order smoothing re- 
sult, Eq. ( |2TI ) is valid, we should expect the ratio of turbulent to 
molecular viscosity to be proportional to the Reynolds number, 

vx^jv Re. (22) 

We find this scaling to be valid for large enough Re in the sim- 
ulation results of Sets A and B, see the upper panel of Fig. |2] 
The vorticity generation is more vigorous for smaller Reynolds 
numbers which can explain the deviating values for small Re. 
The vorticity dynamo exists at least up to Re « 100 (Kapyla 
et al. 120091 ), but the growth rate seems to decrease somewhat as 
Re is increased. We also find that the turbulence becomes more 
isotropic as Re increases (cf. Table IB. ft and conjecture that this 
behaviour is also linked to the less efficient large-scale vorticity 
generation. 

Furthermore, we can estimate the Strouhal number by com- 
paring the simulation results for vx. and the FOSA estimate, 
Eq. ( |2TI ). see the lower panel of Fig. |2] We find that for large 
Reynolds numbers, t't/t'to ~ 1-5, indicating that St w 1.5. 
This value is in accordance with earlier numerical studies of pas- 
sive scalar transport (Brandenburg et al. |20041 l and nondiffusive 
Reynolds stresses (Kapyla & Brandenburg '2008l l. 

The shear dependence of the non-zero components of the 
normalized Reynolds stresses, = Qij/u^^^, from Sets C 
and D is shown in Fig.|3](see also Table IB. il l. We find that Q^x 
decreases and Qyy increases as a function of Sh, although the 
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Sh — 0.2 (Set B) 
Sh — 0.1 (Set A) 



10 



100 



Re 



II 

-l-> 
t/3 




10 



100 



Re 



Fig. 2. Upper panel: turbulent viscosity divided by the molecu- 
lar viscosity for the simulation Sets A (solid line) and B (dashed 
line) as a function of the Reynolds number. The dotted line, 
proportional to Rc, is shown for reference. Lower panel: the 
Strouhal number, St — i/t/t'to, for the same sets of simulations 
as in the upper panel. 

latter trend is barely significant due to the large error bars. Qzz, 
on the other hand, seems to increase slightly for strong shear but 
the error bars are again so large that this trend is not statistically 
significant. For — Sh <; 0.2 the errors are generally quite large 
due to the short averaging interval that we are forced to use due 
to the vorticity generation. 

We can rewrite Eq. iT% in terms of Sh and St by substituting 
vto — ^StUinisfef^^ and using the definition of Sh to obtain 

Q^y = -iStSh. (23) 

This relation is valid for weak shear and is also borne out of 
the closure model (see Appendix IA.lt . We find that the stress 
component Qxy in the Sets C and D is approximately consistent 
with a relation —Qxy oc (0.5 . . . 0.67)Sh, suggesting that St « 
1.5 ... 2 (see the bottom panel of Fig. [3] and the next section). 
The Reynolds number in the Sets C and D varies between 27 
and 123 so the results for the Strouhal number are consistent 
with the corresponding data in Fig. |2] 

3.1.2. Closure model results 

We now turn to the simple closure model that was introduced 
in Sect. 12.41 We compare the stationary solutions of the MTA- 
model with the time and volume averaged simulation data. The 
closure model predicts that the absolute values of Q^x and Qzz 
remain constant as functions of Sh (see Appendix lAl for more 
details) 





0.34 - 


CO 


U.OC 7 


CM 1 


0.30 7 


\ 
»< 


0.28 7 








0.26 7 




0.24 7 



o Re~27 



St=2.5 

St=2 

^ Re~75 st=1.5 



0.00 0.05 0.10 0.15 0.20 0.25 




0.00 0.05 0.10 0.15 0.20 0.25 




0.00 0.05 0.10 0.15 0.20 0.25 




)(0) 

i/ XX ' 



Qzz = Q 



(0) 
zz ■ 



(24) 



0.00 0.05 0.10 0.15 0.20 0.25 
-Sh 

Fig. 3. From top to bottom: Reynolds stress components Qxx, 
Qyy, Qzz, and Qxy, normalised by Ur^js, as functions of the 
shear parameter for simulation Sets C and D with Reynolds num- 
bers, based on Umis from a run with Co = Sh = 0, of ~ 27 (di- 
amonds) and w 75 (triangles), respectively. The curves in each 
panel show the results of the MTA-closure with three Strouhal 
numbers. The linestyles are as indicated in the legend in the top 
panel. 



However, since Qyy, and thus Q — u^^^^, can vary, the nor- 
malised data shows a decreasing trend for Qxx and Qzz as a 
function of shear The MTA-model also shows that Qxy is lin- 
early proportional to the shear parameter for weak shear whereas 
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the other two off-diagonal components are zero as expected from 
symmetry cosiderations. 

We find that the closure model is in qualititive agreement 
with the numerical results for Qxx, Qyy, and Qxy for all val- 
ues of Sh, see Fig. |3] For Q^z the numerical data shows large 
scatter for increasing values of |Sh|, although the points around 
— Sh w 0.12 are consistent with a declining trend. Considering 
the other components, the best fit to the simulation data for Qyy 
and Qxy is obtained if St w 1.5 ... 2 is used in the MTA-model 
(see the dotted and dashed lines in Fig.O, whereas a somewhat 
larger St fits the results of Qxx best. Especially the off-diagonal 
component Qxy is quite well reproduced with our simple model. 

3.2. Case n ^ 0, 9 = 

When adding shear, we consider two distinct lines in the param- 
eter space: in simulation Set E we keep Sh constant and vary 
Co, whereas in Set F the rotation is kept constant and the shear 
is varied (see also Table iBTZt . The range in which these param- 
eters is varied is given by —5 < q < 1.9 where the parameter 
q = —S/flo can be considered to describe the shear in a differ- 
entially rotating disk (see Sect. 12.3b . Whilst the cases q — 1.5 
and q = 1 can be thought to represent Keplerian and galactic 
disks, respectively, the regime q < 1 is not likely to occur in the 
bulk of the disk in any system, but such configurations can oc- 
cur in convectively unstable parts of stellar interiors, such as the 
solar convection zone. Cases g > 2 are Rayleigh unstable (see 
Appendix IB. lb and lead to a large-scale instability of the shear 
flow so simulations in this parameter regime are not considered 
here. 



3.2.1. Simulation results 

In Set E we fix the shear flow, measured by Sh, and vary the ro- 
tation rate. In this case the value of Sh w —0.16 is very close to 
being constant in all runs. The simulation results for the nonzero 
components of the stress are shown in Fig. |4] We find that the 
components Qxx and Qyy behave very similarly to each other 
but with opposite trends as functions of q, i.e. the sum of the two 
is roughly constant. The quantity Q^z varies much less than the 
two other diagonal components, and is consistent with a constant 
value as a function of q within the error bars in all cases apart 
from three points in the range \q\ < 0.5. The off-diagonal stress 
Qxy shows almost a symmetric profile with respect to q = 0, 
with somewhat steeper rise on the positive side. We note that 
since Sh w const < 0, the very small values of Qxy around 
q = indicate that the more rapid rotation there either quenches 
the turbulent viscosity or that additional nondiffusive contribu- 
tions with the opposite sign are present. 

In Set F, we keep constant and vary the magnitude of 
the shear In this case, however, the Coriolis number no longer 
stays exactly constant because the rms-velocity is affected by the 
rather strong shear for the extreme values of\q\. This also affects 
the anisotropy of the turbulence which is much greater than in 
Set E, see Fig.|5] The trends of Qxx and Qyy are again opposite 
to each other as in Set E, but here the variation as a function of q 
is significantly greater Similarly as in Set E, the Q^^-component 
is less affected, although a weakly incresing trend is seen for 
small q and a clearly decreasing trend is seen for positive values 
of q. The off-diagonal component Qxy changes sign at q — Q 
where also Sh changes sign. However, it is clear that Qxy is not 
linearly proportional to q as might naively be expected if the 
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Fig. 4. Same as Fig. |3]but for Set E, and different Strouhal num- 
bers as indicated by the legend in the uppermost panel. 



stress is fully due to turbulent viscosity. Thus it is important to 
try to sepatate the nondiffusive and diffusive contributions (see 
below). 

3.2.2. Closure model results 

The Reynolds stresses obtained from the closure model are com- 
pared with the results of the numerical simulations from Sets E 
and F in Figs. |4] and |5] with constant Sh and Co, respectively. 
In both cases we observe rather good qualitative agreement be- 
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Fig. 6. Upper panel: the two measures of the horizontal A-effect 
from Set E according to Eq. (|29] l (solid line) and Eq. ( [32] i 
(dashed). Lower panel: the same as above but from Set F. 



studied here. In Set E the data is within error bars consistent 
with St = 0.75 . . . 1.00 in the range q < 0.5 except for Qzz near 
q = 0. Furthermore, St > 1 would be required to reproduce Q^x 
and Qyy for 0.5 < g < 2. In Set F the situation is similar with 
Qzz and Q^y deviating for large negative q and Qyy deviating 
for g > 1. All in all, the correspondence between the simple 
closure and the simulation results is rather good. 
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Fig. 5. Same as Fig.|4]but for the simulation Set F. 



tween the numerical simulations and the closure model with the 
exception of Qzz in Set F when q < —2 (see the third panel 
of Fig. |5j. Another discrepancy for Qzz occurs in the Set E 
near q = where the simulation results indicate a higher value 
than what is obtained from the closure model. The discrepancies 
for Qzz occur for rapid rotation (Set E near q — Q) and large 
shear (Set F for q < — 2) which suggests that the simple closure 
model used here could be inadequate in such regimes. We also 
note that the fit generally tends to get worse as q approaches 2. 
Partially due to this we cannot ascribe a single Strouhal num- 
ber that would fit the simulation data in the full parameter range 



3.2.3. Separating diffusive and nondiffusive contributions 

When both, shear and rotation are present, the Reynolds stress 
can also contain nondiffusive contributions proportional to the 
rotation rate, i.e. 

Qij — Mjk^k — J^ijklUkd, (25) 

where Aijk describes the A-effect (Krause & Riidiger 119741 
Rudiger 1980l. [T989l Kitchatinov & Rudiger [T993] l, which con- 
tributes to the generation of differential rotation. In the case that 
we have shear in addition to rotation, the problem is that for 
a given stress component we have two unknown coefficients 
which makes it difficult to distinguish the different contributions 
without a method similar to the test field procedure in magneto- 
hydrodynamics (see, e.g. Schrinner et al. ' 20051120071 ). However, 
if we consider that the diffusive contribution is given by Eq. ( fT9] l 
with i^t ~ St h'to^ it is possible to separate the two provided that 
the shear and rotation are sufficiently weak. 

Considering the case ft — rio(0, 0, 1), the nondiffusive part 
of the Reynolds stress can be described by a single coefficient 
that is commonly denoted (cf. Rudiger[l989|l by 

gi^) = Ancos^rJo. (26) 

This is often referred to as the horizontal A-effect because it 
generates horizontal differential rotation in mean-field models 
of stellar interior rotation. Let us now write the total stress as 
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a sum of the contributions from the A-effect and the turbulent 
viscosity 

Qxy = Anr^o - ^tS. (27) 

The turbulent viscosity is given by vt — ^tu^^-^^ with r = 

St/(ui.nisfcf), so = iSt-Urms^f^^ = Stfc'to- Substituting this 
above we obtain 

Qxy = Ah^^O " StftoS*. (28) 

Solving for Ah and dividing by lyto^o yields 
AH^AW=6|^~Stg, (29) 

where Q^y = Qxy/u^rns^ Ah = An/i^to, q = -S'/fio and 
where the definition of Co, Eq. (|9]l has been used. 

On the other hand, for slow rotation and in the absence of 
shear the nondiffusive stress due to the anisotropy of the turbu- 
lence can be written as (Riidiger l 19891 see also KB08) 

Qi^' =217,t(Q^^-Q,,). (30) 
Equating this with Eq. (|26] | gives 

Kll=2T{Qyy-Q^^). (31) 

Dividing by i/to and eliminating r using St yields 

Ah = A^j') = 6 St{Qyy - Q,,). (32) 

We can now compare the expressions ( |29] l and (|32] | by using the 
stresses and Co from the simulations and keeping the Strouhal 
number as a free parameter. 

The results for the two measures of the A-effect for the Sets E 
and F are shown in Fig. |6] In Set E the qualitative behaviour of 
the two expressions is the same for Strouhal numbers greater 
than « 0.7, whereas the best correspondence between A^^"* and 
Ay is obtained for St w 0.85. For larger St the qualitative trend 

of Ay^ stays the same but the deviation between the two ex- 
pressions becomes increasingly greater for large values of In 
summary, in Set E we seem to be able to extract a reasonable 
estimate of the nondiffusive contribution to the Reynolds stress 
with the method outined above. The magnitude of Ah is of the 
order of (0.2 . . . 0.4)//to when a similar Strouhal number is used 
in the fitting as was required for the closure model to approxi- 
mately reproduce the simulation results. The magnitude of Ah 
is also quite close to the values obtained by KB08 in a system 
without shear. 

In Set F the qualitative behaviours of and A^^"* are the 
same for Strouhal numbers greater than about w 0.5. We find 
that in the Set F a g-independent Strouhal number does not give 
a very good fit to the numerical data. We have plotted the curves 
for St = 0.75 which yields a reasonable fit to the data in the 
vicinity of g = where |Sh| is small, and Eqs. ( |29] l and ( [32] l 
can be considered to be the least affected by the shear. Thus the 
results for the A-effect in Set F leave more room for speculation. 
One contributing factor is likely to be the significantly stronger 
shear for the extreme values of \q\ which possibly renders the 
expressions ( |29l ) and ([32]) inaccurate in this regime. 



3.3. Case 0,9 = 90° 

Finally, we consider the case where, in addition to the shear flow, 
rotation with 6 = 90° is used, i.e. = Jlo(— 1, 0.0). As is noted 
by Leprovost & Kim (2 008 b) this system is unstable even in 
the absence of turbulence. The instability generates large-scale 
vorticity due to the fact that no stationary solution to the equa- 
tions of the large-scale velocity exist (Leprovost & Kim 2008bl 
see also Appendix |B.2t in the absence of forcing on the large 
scales. If the system is interpreted as the polar region in the so- 
lar tachocline, such large-scale forcing can be available, e.g. by 
thermally driven flows but in the present simulation setup such 
effects are difficult to implement. 

In spite of the unstable nature of the system, we find that 
when turbulence due to the forcing is present we can still extract 
information about the Reynolds stresses from the early stages of 
the simulations where the large-scale flows are still weak. The 
procedure is similar to the case where rotation was absent (see 
Sect. I3.1.U . However, this means that the error bars in Figs. Q 
and |8] are greater than those in Figs. |4] and |5] because much 
shorter averaging intervals must be used. The situation is ag- 
gravated in Set H (see, Fig.fHJ where the amplitude of the shear 
flow increases when \q\ increases. 

3.3.1. Simulation results 

We perform again two sets of simulations (see Table IB. 3"T l where 
either the shear flow (Set G) or rotation (Set H) is kept fixed. In 
the previous cases the source of anisotropy was either only large- 
scale shear with mean vorticity W = V x J7 = Sz (Sets A to 
D), or shear and rotation of the form ft = flQZ (Sets E and 
F). In the present case the rotation vector points to the negative 
x-direction due to which further symmetries are broken and the 
other off-diagonal stresses Q^z and Qyz can have nonzero values 
(e.g. Rudiger 1989). 

The results from Set G are shown in Fig.|2] Due to the large 
error bars, the results for the diagonal components are consis- 
tent with a value independent of q. The component Q^x seems 
to increase and Qyy to decrease near q = 0, but these results are 
not statistically significant. The off-diagonal stress Q^y shows a 
similar qualitative and quantitative trend as in Set E. The compo- 
nent Qxz is negative (positive) for q < (q > 0) with absolute 
values peaking near q = 0. The magnitude of Q^z is approxi- 
mately one third of Qxy. The Qyz component is smaller by an- 
other factor of two to three in comparison to Q^z with a profile 
consistent with linear proportionality to q. However, the large er- 
rors for q < 0.25 imply that the results there are not statistically 
significant. 

The results from Set H are shown in Fig. [8] The diagonal 
stresses show again little dependence on q except for q < —2.5, 
where Q^x decreases and Qyy increases, although the error bars 
are again so large that the results for Qyy are also consistent 
with a constant profile. Within error estimates, Qzz is consistent 
with a constant profile as a function of q. In general, the results 
for q < —2.5 are rather unreliable due to the rapid generation 
of large-scale flows. The off-diagonal stress Q^y has a similar 
profile and somewhat larger magnitude than in Set F. Since the 
shear flow in both sets is the same the difference can be due to 
the different nondiffusive contributions. Component Q^z shows 
a similar profile, and a magnitude of about one third of the stress 
Qxy Qyz, on the other hand, exhibits a profile symmetric around 
q = 0, apart from a few unrealiable points in the regime q < 
-2.5. 
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Fig. 7. Left column: Qxx (top panel), Qyy (middle), and Qzz (bottom) from the simulation Set G. Right column: Qxy (top panel), 
Qxz (middle), and Qyz (bottom) from Set G. The curves in each panel show the closure model results with three different Strouhal 
numbers as indicated in the legend (righ column, middle panel). 



3.3.2. Closure model results 

The Reynolds stresses obtained from the closure model are com- 
pared with the results from the numerical simulations of the 
Sets G and H in Figs. [7] and [8] respectively. In Set G the error 
bars for the diagonal components are so large that a wide range 
of Strouhal numbers are consistent with the results. Because the 
shear in this set is rather weak in all runs, the anisotropy of the 
turbulence remains weak in all simulations. It also appers that 
the simulation results for Qxx and Qzz are not captured by the 
closure model for some points near q ~ 0. The off-diagonal 
components Qxy and Qxz show a clearer picture and the closure 
model is consistent with the numerical results if St = 0.75 ... 1. 
The simulation results for Qyz are for the most part not distin- 
guishable from zero for q < 0.5. For q greater than this, the 
closure model is consistent with the simulations for St = 1. 

In Set H the diagonal components show weak signs of 
anisotropy if \q\ < 2. However, the errors increase substantially 
for greater q and the data is still almost consistent with isotropic 



turbulence even for q = —5. This is due to the short averaging 
interval that results in from the vigorous vorticity generation that 
is excited very early in these runs. This is because the extreme 
values of \q\ correspond to large values of Sh in this set. All of 
the off-diagonal stresses show qualitatively correct behaviours 
and the simulation data is again consistent with St = 0.75 ... 1 
for |g| < 3. 

3.3.3. Separating diffusive and nondiffusive contributions 

In the case with 6 = 90°, Eq. ( |26] l no longer applies, but we can 

still use Eq. (|29] l to estimate the nondiffusive contribution A[]^\ 
On the other hand Eq. ( |32] | is no longer valid because flz = 0. To 
replace Eq. ( [32] i, we can use the tau-approximation to derive an- 
other expression for the nondiffusive stress. We arrive at a term 
proportional to fix in the equation of Qxy 

Q<^1^ =2nxTQxz = A'i^no, (33) 
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where the superscript differentiates this contribution from those 
used in the case where 9 ~ Q. Division by i^tQ yields 



A 



(3) 



(34) 



Furthermore, to lowest order, the other two off-diagonal stresses 
can contain nondiffusive contributions 

Qi^' = Am^o, (35) 
Q^^) = AvsinSf^o- (36) 
Solving for the A-coefficients and normalizing by i^to gives 



from symmetry considerations in a rotating spherical system. 
However, here these symmetry considerations are no longer 
obeyed due to the added shear flow. Furthermore, applying 
the minimal tau-approximation to the Reynolds stress equation 
yields 



2 ^x^Q xy 1 



'vv 



)• 



(39) 
(40) 



Equating these with Eqs. (|35] |. and ( |36] |. respectively, and bear- 
ing in mind that 0,^ — — 57o, we obtain, after normalization with 
the turbulent viscosity vto. 



Am = Aj^ — 



Av ^ A« = 



Co ' 
Co ■ 



(37) 



(38) 



A^^ — SStQajy, 



(41) 
(42) 



We note that whilst Eq. ( [36] l coincides with that commonly 
assumed for the A-effect, Eq. ( [35l l does not. Typically the 
"meridional" A-effect, Am, is proportional to sin 6* cos 6*51 (e.g. 
Pulkkinen et al. 119931 Kapyla et al. i2004i which derives 



We now proceed with the comparison in the same manner as in 
Sect. 1123] 

We find that in Set G the correspondence between the two 
measures for each of the A-coefficients are in rough qualitative 
agreement for St w 1, but matching the quantitative values is not 
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Fig. 9. Nondiffusive contributions to the Reynolds stresses pa- 
rameterized by the coefficients Ah (top panel), Am (middle), 
and Ay (bottom) from the simulation Set G. In each panel two 
expressions are compared and the linestyles are explained in the 
legends. 

possible (see Fig.|9]l. Compared to Set E, the profile and ampli- 
tude of Ah ~ 0.05 . . . 0.3 are similar, although a precise value 
is difficult to determine due to the disagreement between the 
two expressions. The other two nondiffusive components have 
magnitudes similar to Ah at least in the range \q\ < 2. For Am 
the analytical expressions seem to agree rather well in the range 
l^l < 2. For Ay, on the other hand, the correspondence is not 
very good, although the sign of the coefficient is for the most 
part reproduced correctly. Here, however, the indeterminate na- 
ture of Qyz is likely to affect the results from Eq. ( [38T l. 

In Set H the correspondence between the different expres- 
sions for the A-coefficients is somewhat better. For g < — 2 the 
simulations begin to exhibit vigorous vorticity generation at a 
very early stage making the error estimates increase significantly 
and the results for the A-effect in this regime are not very re- 
Uable. Again the best fit is achieved with St w 1, although a 
unique value that would fit all components cannot be found. The 
maxima of all A-coefficients are of the order 0.5 in this set. 

3.4. On the validity of tlie MTA 

We note that the largest deviations between the simulations and 
the closure model are seen when rapid rotation or strong shear 
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are used. This may suggest that the simple closure model used 
here can break down at such circumstances. In order to test the 
validity of the present formulation of the MTA, we perform a 
set of simulations where q ~ 1 is constant and the values of 
Co and Sh are varied. These results are compared to those from 
the MTA-closure model where the relaxation time is indepen- 
dent of Co and Sh. The results are shown in Fig. [TT] It appears 
that the simulation results for all of the stress components are 
in accordance with the closure model only for Co ^ 0.1. The 
components Q^x and Qyy are not well reproduced for greater 
Co. Rather surprisingly the Q^z component is consistent with 
the closure model for all values of Co explored here. For Q^y 
the agreement is rather good for Co ^ 0.3, but for greater Co 
the stress changes sign as a function of Co in the simulations 
which is not reproduced by the closure model. Due to the dif- 
fering behaviour of the stresses we cannot estimate the validity 
range of the MTA-closure very precisely based on this data only. 
It is likely that with strong enough rotation or shear the relevant 
time scale to be used in the closure model is associated with ro- 
tation or shear instead of our naive estimate of the eddy turnover 
time Eq. ( fT6b . For comparison, the time scales related to rotation 
and shear are of the same order of magnitude as the turnover time 
of the turbulence when Co — 0.5. At face value, the present re- 
sults suggest that the effects of rotation and shear begin to affect 
the results even at somewhat slower rotation. 

Furthermore, it turns out that in the presence of a shear flow 
of the form Uy{x) it is possible to estimate the relaxation time 
using the basic assumption of the MTA, Eq. ( fT3] ), from the equa- 
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tion of Qxy This is analogous to the passive scalar case where 
the flux due to an imposed gradient was found to the propor- 
tional to the triple correlation (Brandenburg et al. 120041) . Let us 
now write 



Qxy 
J^xy 



where 



Txy — —CgUidj Inp + Ujdi Inp + QijV ■ u. 



(43) 



(44) 
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Fig. 12. Upper panel: Strouhal number from the triple correla- 
tions as a function of Co for the same runs as in Fig. ( fTTT i. Lower 
panel: St as a function of Re for runs with Co ~ 0.06 and q = 1. 

We find that the last term of Eq. ( |44] | is negligible in comparison 
to the other terms in all cases considered here. Representative re- 
sults for the Strouhal number defined by Eq. ( fT6] l from two sets 
of runs are shown in Fig. [12] Firstly, as a function of rotation, the 
Strouhal number is between 1.5 and 3 for Co ^ 0.3, whereas the 
value decreases for more rapid rotation and shear. The relaxation 
time can even be negative for Co ^ 1, indicating that the basic 
assumption of the MTA breaks down. Furthermore, for weak ro- 
tation and shear, the Strouhal number is constant as a function of 
the Reynolds number provided that Re > 1. This is consistent 
with the fact that in the regime Re w 1 the viscous time scale is 
of the same order of magnitude as the turnover time and neglect- 
ing the viscous terms in the equation of the stress is not justified. 
However, as we use Reynolds numbers of the order of 30 in most 
of our simulations, omitting the viscosity is permitted. 

4. Conclusions 

We have performed simulations of isotropically forced homoge- 
neous turbulence under the influences of shear and rotation in 
order to study the turbulent transport properties. We find that in 
the absence of rotation the turbulent viscosity is of the order of 
the first order smoothing estimate i^to = \u^-mskY^, and that 
the Strouhal number saturates to a value of St w 1.5 for large 
Reynolds numbers. 

When both rotation and shear are present, we seek to distin- 
guish the diffusive (turbulent viscosity) and nondiffusive parts 
(A-effect) of the stress. In the case of Q^y this is achieved by 
assuming the turbulent viscosity to be proportional to Stz/to and 
solving for the A-part from the equation of total stress, whereas 
for the other two off-diagonal components no diffusive contri- 
butions are assumed to appear in the equation of the total stress. 
As an independent check we derive analytical expressions for 
the A-effect using the minimal tau-approximation. The two ex- 
pressions are then compared by using the Strouhal number as 
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a free parameter. Although the results for the A-effect from the 
different approaches are not fully consistent in all cases, we find 
that they agree at least qualitatively, if not quantitatively, in most 
cases. The magnitude of the A-effect is a few times O.lz^to in 
most cases in accordance with earlier studies (e.g. Kapyla & 
Brandenburg 120071 120081) . However, a more precise determina- 
tion and separation of the different coefficients requires a method 
akin to the test field method used succesfully in magnetohydro- 
dynamics (Schrinner et al . 120051120071 1. 

We have also studied the Reynolds stresses using the mini- 
mal tau-approximation closure model where the Strouhal num- 
ber St is the only free parameter Comparing the results of this 
very simple closure model with our simulations shows that they 
generally agree surprisingly well with rather few exceptions. We 
find that in most cases the best fit between the closure model 
and full simulations is achieved when St w 1. We note, how- 
ever, that the agreement between the closure model and the sim- 
ulations becomes worse as the rotation rotation rate or shear are 
increased. The reason is likely to be that the relevant time scale 
in those regimes is no longer the turnover time but rather the time 
scale related to the rotation or shear. We feel that investigating 
the validity of the tau-approximation more precisely and futher 
improving its performance in reproducing the simulation results 
is a valid subject for further study. 

Finding a closure model that reproduces the relevant prop- 
erties of turbulence would be very useful in a variety of astro- 
physical applications. From the viewpoint of the present paper 
the most interesting application is the angular momentum bal- 
ance of convectively unstable stellar interiors where rotation and 
shear flows are known to play important roles. Although the 
present model is highly idealised in comparison to stratified con- 
vection, the present results using the MTA-closure are promis- 
ing. A logical step towards more general description of con- 
vective turbulent transport would be to study Boussinesq con- 
vection (e.g. Spiegel & Veronis 1960) and extend the closure 
model to work in the same regime (Miller & Garaud 120071 ). 
In addition to the Reynolds stresses, such models could be 
used to study anisotropic turbulent heat transport (e.g. Riidiger 
[19821 Kitchatinov et al. [19941 Kleeorin & Rogachevskii 2006). 
Furthermore, these models could help in answering how some 
enigmatic results of local convection simulations, such as the 
narrow regions of very strong stresses near the equator (Chan 
[2001 ; Kapyla et al. [2004l Rudiger et al. '2005') in the rapid rota- 
tion regime, can be explained. These subjects will be studied in 
subsequent publications. 
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Appendix A: Stationary soiutions to tlie 
MTA-closure modei 

Stationary solution to Eq. dTSl ) can be obtained by setting Qij — 
and solving the resulting equations 



^xx — ^^z^Qxy ~t~ Qxx ^ 

2xy StQxx 2 ^x^Q xz 2 ^z^(,Qyy 

^xz — ^ ^z'^Qyz ^^x^Qxy-) 



Qxx\ 



(A.l) 
(A.2) 
(A.3) 



Qyy = -2 SrQ^y + 4 n^rQy,^ - 4 QzTQxy + Q'y^J , (A-4) 
Qyz =- -{S + 2 n,)TQ,, + 2 n,T{Q,, - Qyy), (A.5) 
Q,, = -4.n,TQy,+Q<-°], (A.6) 
where ft^ and ft^ are the components of the rotation vector 
as defined in Eq. (ITOl ). and we have used the fact that Qi""* = 

Qxz — Qyz — 0. 



A. 1. Only shear, Sh 7^ 0, Co = 0. 

In the absence of rotation, the stationary solution to Eqs. ( lA.l 



to (IA.6b 


is 


Qxx — 


^xx ) 


Qxy 


-StShQ( 


Qyy — 


2St2Sh2( 


Qzz = 


^ZZ 1 


Qxz — 


Qyz — O7 


Q - 





(0) 

XX "! 



^yy ■ 



(A.7) 
(A.8) 
(A.9) 

(A. 10) 
(A. 11) 
(A. 12) 

Since the forcing in the present case is isotropic, we can write 

Qxx Qyy Qzz sQ^ ^' (A. 13) 

where Q^"' is the trace of Q^^'' in the case without rotation or 
shear Substituting (IA.13I ) into the expression of Q, we obtain 



]^ + Q^- 



= Q = Q'"^ (1 + ISt^Sh^ 



(A. 14) 



Using Eq. (|A.13t and ( IA.14I ) the normalised stresses, Qij = 
Qij/u^ais '^^n be written in terms of the dimensionless param- 
eters St and Sh as 



Qxx — Q 

Qxy 



1 



3 + 2St2Sh2^ 
-StSh 



3 + 2St2Sh2' 



■ yy 



= 1 - 



3 + 2St2Sh^' 



(A.15) 
(A. 16) 
(A.17) 



A.2. Shear and rotation, Sh, Co 7^ 0, 6i = 0. 

If rotation is present in the system and 9 = 0, the solution of 
Eqs. (lA.lb to (IA.6b can be written as 

[2CoSt2(Co + Sh) + l]Qi*L^+2Co2St2g|,°) 

Q^X — : . ^ ?/^ ' (A.l«) 



Q 



l + 4CoSt2(Co + Sh) 
-(Co + Sh)StQi"i +CoStQ^ 



)(0) 



xy 



l + 4CoSt^(Co + Sh) 



(A.19) 



[2 CoSt^(Co + Sh) + l]Q 
l + 4CoSt^(Co + Sh) 
2 (Co + Sh)2Qi°j 



(0) 



l + 4CoSt2(Co + Sh)^ 
1(0) 



Qxz — Qyz — 0; 

Q = (i 



2Sh2St2 



Q 



(0) 



(A.20) 

(A.21) 
(A.22) 

(A.23) 



3 + 12CoSt2(Co + Sh), 

If we allow r 7^ Tf the solution can be obtained by simply mul- 
tiplying the above solutions with a factor r/rf. 
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A3. Shear and rotation, Sh, Co^O,e = n/2. 
Considering the case 9 = tt/2, the solution is 

Qxx^ 



ixx — ^xx 



Qxy 



Qyy — 



Ql°ishSt 
1 + Co^St^ ' 
Qi"i St" CoSh 
l + Co"St" ' 
2 Sh'St'Qi°i+Q<"J +2 (g|,°J + Qi"i)Co'St' 
l + 4Co"St" 
-Qi°i)Co'St']CoSt 
(l + 4Co2St2)(l + Co2St2) 



(A.24) 
(A.25) 

(A.26) 

,(A.27) 



(l + 4Co'St')(l + Co^St") 

2 (Q^°,^ + Q£^)Co^St^ + 6St^gi°iCo^Sh^ 
(l + 4Co2St2)(l + Co^St2) 



(A.28) 



(l + 4Co^St^)(l + Co^St^) 



Q ^ 1+ o 



2 Sh^St^ 



3 1 + Co^St^ 



(A.29) 
(A.30) 



If T 7^ T/, the solution can be obtained in the same way as in the 
= case. 



Appendix B: Stability analysis 

B. 1. Vorticity dynamo in tlie presence of rotation 

The existence of the vorticity dynamo in the present case can be 
studied using the same procedure used in Kapyla et al. (2009). 
When one assumes that the average velocity field depends only 
on z, the relevant mean field equation has the form 



U = -SUxV -u-Vu-2nxU + i'U , 



(B.l) 



where the primes denote z-derivatives. Assuming the flow to be 
incompressible, V • J7 = Uz.z = 0, so Uz = const = by 
a suitable choice of the initial condition. The nonlinear term 
—u ■ Vm can be represented in terms of the large-scale veloc- 
ity by introducing a turbulent eddy viscosity tensor u^j (see, e.g. 
Elperin et al. 120031 Kapyla et al. l2009l l 



- {u ■ Vu)i = UijUj . 
Substituting Eq. ( IB.2b to Eq. (IB. Il l gives 

Ux ^ {U + Vxx)Ux" + VxylJy" + 2^lzUy, 

Uy = [U + Vyy)Uy" + UyxUx" - {S + 2Vlz)Ux- 



(B.2) 



(B.3) 
(B.4) 



We note that because U z — 0, no terms proportional to 0,^ ap- 
pear Similarly as in Kapyla et al. 2009 , we assume a stationary 
state and seek solutions to Eqs. (IB. 3b and (IB.4b . Now we can 
apply a Fourier transform and arrive at 



{v + Vxx)k^Ux + {vxyk^ - 2 nz)(jy = 0, 

{V + Vyy)k^Uy + [Vyxk^ +S + 2 Vlz^x = 0, 



(B.5) 
(B.6) 



where Ui are the Fourier amplitudes of and k is the 
wavenumber Equations ( IB.5I ) and ( IB.6I ) imply that the sufficient 
condition for the mean vorticity dynamo to appear is 



{l^xy 2-p-)(l^ya; + 



> 1, 



(B.7) 



where vt = ^ + ^{ly^x + >^yy) and e = ^{v^x - Vyy). In the 
absence of shear we recover the expression derived in Kapyla et 
al. 120091 suggesting that v^yS > 0. On the other hand, if we set 
u = and i^ij = 0, we have 

2^(<Z-2)>0, (B.8) 

where (7 = —S/ilz-Since fi^/fc^ is always positive, this inequal- 
ity holds if q > 2. This is the standard Rayleigh stability con- 
dition (for instability). It is clear that rotation has a stabilizing 
effect in the system provided that g < 2. It is, however, con- 
ceivable that for small enough q the component j^xy overcomes 
the rotational stabilization. However, even for our most extreme 
runs with q = —5 this appears not to be the case. 

B.2. Force balance for Sh, Co ^ 0, 6 = 7r/2 

We find that large-scale flows tend to always appear in the 
simulations at = tt/2. A very similar setup was studied by 
Leprovost & Kim ( 1200 8b ) who studied the force balance for the 
large-scale flow including also the pressure gradient 

U = -SUxV -u-Vu-2nxU --VP + vTj" . (B.9) 

P 



Assuming that there is no turbulence (— m • Vm = 0) and no 

large-scale flows {U = 0), apart from the imposed shear U^"^ = 
{0, Sx,Q), one finds the following equations 



dxP = 0, 
dzP-2pnoSx = 0. 



(B.IO) 
(B.ll) 



We see that Eq. dB.lOb indicates that the pressure P is indepen- 
dent of X. The other equation dB.llb . however, gives an explicit 
x-dependence for P. Therefore, we can conclude that unless 
there exists some sort of large-scale forcing to cancel the pres- 
sure gradient, there can be no equilibrium for this system. As 
noted by Leprovost & Kim ( 2008b ), in the real physical situa- 
tions phenomena like the thermal winds can balance the system. 
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Table B.l. Summary of the runs with shear flow only. Here Qij = Qij /u^^^, and brackets around a quantity signify that the number 
is not statisticaUy significant. 



Run 


Grid 


Re 


Co 


Sh 


Urms 


Qxx 


Qyy 


Qzz 


Qxy 


Qxz 




Qyz 


Al 


64^ 


5 





-0.18 


0.054 


0.295 


0.404 


0.302 


0.075 


(-0.000) 


( 


-0.001) 


A2 


128^ 


14 





-0.14 


0.069 


0.304 


0.378 


0.319 


0.084 


(-0.001) 




(0.002) 


A3 


256^ 


31 





-0.12 


0.079 


0.313 


0.368 


0.320 


0.083 


(-0.002) 


( 


-0.001) 


A4 


256^ 


80 





-0.12 


0.083 


0.331 


0.342 


0.328 


0.069 


(-0.001) 


( 


-0.001) 


A5 


512^ 


164 





-0.12 


0.084 


0.318 


0.349 


0.333 


0.058 


(-0.003) 


( 


-0.003) 


A6 


1024^ 


345 





-0.11 


0.088 


0.320 


0.343 


0.337 


0.054 


(-0.001) 


( 


-0.003) 



Bl 


64'^ 


8 





-0.25 


0.080 


0.198 


0.599 


0.238 


0.146 


(0.002) 




(0.001) 


B2 


128'^ 


22 





-0.18 


0.112 


0.235 


0.491 


0.288 


0.140 ( 


-0.001) 


( 


-0.000) 


B3 


256'^ 


44 





-0.18 


0.112 


0.262 


0.440 


0.303 


0.124 ( 


-0.002) 




(0.000) 


B4 


512^ 


99 





-0.19 


0.101 


0.282 


0.425 


0.318 


0.110 


(0.001) 




(0.001) 


B5 


1024^ 


196 





-0.20 


0.100 


0.300 


0.377 


0.323 


0.101 ( 


-0.003) 


( 


-0.001) 


B6 


1024^ 


387 





-0.20 


0.099 


0.299 


0.375 


0.326 


0.095 ( 


-0.004) 


( 


-0.005) 


CI 


128^ 


27 





-0.01 


0.069 


0.335 


0.334 


0.332 


0.010 ( 


-0.003) 


( 


-0.002) 


C2 


128^ 


27 





-0.03 


0.069 


0.334 


0.332 


0.335 


0.014 ( 


-0.001) 




(0.001) 


C3 


128^ 


28 





-0.06 


0.070 


0.335 


0.332 


0.333 


0.031 ( 


-0.002) 




(0.000) 


C4 


128^ 


31 





-0.13 


0.078 


0.317 


0.362 


0.321 


0.084 ( 


-0.001) 


( 


-0.000) 


C5 


128^ 


39 





-0.20 


0.099 


0.282 


0.392 


0.330 


0.125 


(0.002) 


( 


-0.002) 


C6 


128^ 


48 





-0.24 


0.122 


0.259 


0.396 


0.347 


0.137 


(0.000) 




(0.010) 


Dl 


256^ 


75 





-0.01 


0.076 


0.331 


0.334 


0.335 


0.009 


(0.002) 




(0.002) 


D2 


256'^ 


75 





-0.03 


0.076 


0.335 


0.329 


0.336 


0.015 


(0.001) 




(0.000) 


D3 


256'^ 


76 





-0.05 0.077 0.329 0.334 0.337 0.029 ( 


-0.001) 


( 


-0.000) 


D4 


256^ 


81 





-0.12 


0.083 


0.331 


0.342 


0.328 


0.069 


(0.001) 


( 


-0.000) 


D5 


256^ 


99 





-0.19 0.101 0.297 0.373 0.331 0.108 ( 


-0.000) 




(0.001) 


D5 


256^ 


123 





-0.23 


0.126 


0.267 


0.417 


0.317 


0.129 


(0.002) 




(0.010) 
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Table B.2. Summary of the runs in the Sets E and F with ^ = 0. 



Run 


Grid 


Rc 




y 




Co 




Sh 


1 


rms 


Qxx 


Qyy 


Qzz 


Qxy 




^ xz 






El 


128'^ 


24 


-5 


00 


-0 


06 


-0 


16 





062 





324 





348 





328 





043 


(0 


000) 


(-0 


000) 


E2 


128^ 


24 


-4 


00 


-0 


08 


-0 


16 





062 





323 





348 





330 





043 


(0 


000) 


(0 


000) 


E3 


128^ 


24 


-3 


00 


-0 


11 


-0 


16 





062 





321 





350 





329 





041 


(0 


001) 


(0 


001) 


E4 


128^ 


24 


-2 


50 


-0 


13 


-0 


16 





061 





320 





350 





330 





039 


(0 


000) 


(0 


000) 


E5 


128^ 


24 


-2 


00 


-0 


16 


-0 


16 





061 





318 





352 





330 





036 


(0 


001) 


(-0 


000) 


E6 


128^ 


24 


-1 


50 


-0 


21 


-0 


16 





061 





317 





353 





330 





033 


(0 


000) 


(-0 


001) 


E7 


128'^ 


24 


-1 


00 


-0 


32 


-0 


16 





061 





315 





353 





332 





027 


(-0 


000) 


(-0 


000) 


E8 


128^ 


24 


-0 


50 


-0 


64 


-0 


16 





061 





315 





351 





335 





014 


(0 


001) 


(-0 


001) 


E9 


128^ 


24 


-0 


25 


-1 


27 


-0 


16 





062 





318 





341 





341 





005 


(0 


001) 


(-0 


001) 


ElO 


128-'* 


24 





25 


1 


28 


-0 


16 





061 





352 





310 





337 





003 


(-0 


001) 


(-0 


000) 


Ell 


128'^ 


24 





50 





65 


-0 


16 





061 





361 





308 





331 





015 


(-0 


001) 


(-0 


001) 


E12 


128^ 


24 





75 





43 


-0 


16 





061 





361 





309 





330 





026 


(-0 


001) 


(-0 


000) 


E13 


128^ 


24 


1 


00 





32 


— 


16 





061 





360 





311 





329 





032 


(— 


001) 


(— 


001) 


E14 


128^ 


24 


1 


25 





26 


— 


16 





061 





358 





313 





329 





037 


(— 


001) 


(— 


000) 


F15 


128^ 


24 


1 


50 





21 


— 


16 





061 





357 





315 





328 





040 


(0 


000) 


(-0 


001) 


E16 


128^ 


24 


1 


75 





18 


— 


16 





062 





355 





317 





328 





043 


(-0 


001) 


(-0 


001) 


E17 


128^ 


24 


1 


90 





17 


-0 


16 





062 





355 





318 





327 





044 


(-0 


001) 


(-0 


001) 


Fl 


128^ 


30 


-5 


00 





25 





64 





077 





230 





431 





340 


-0 


075 


(-0 


001) 


(0 


000) 


F2 


128^ 


28 


-4 


00 





27 





55 





072 





253 





407 





340 


-0 


068 


(-0 


001) 


(0 


001) 


F3 


128-'* 


26 


-3 


00 





29 





44 





067 





276 





388 





336 


-0 


058 


(-0 


001) 


(0 


002) 


F4 


128'^ 


26 


-2 


50 





30 





37 





066 





287 





381 





333 


-0 


053 


(-0 


001) 


(0 


002) 


F5 


128'^ 


25 


-2 


00 





31 





31 





064 





297 





370 





333 


-0 


046 


(-0 


002) 


(0 


002) 


F6 


128^ 


25 


-1 


50 





31 





24 





063 





306 





363 





330 


-0 


038 


(-0 


001) 


(0 


003) 


F7 


128^ 


24 


-1 


00 





32 





16 





061 





314 





355 





331 


-0 


027 


(-0 


002) 


(0 


004) 


F8 


128^ 


24 


-0 


50 





33 





08 





060 





325 





344 





331 


-0 


013 


(-0 


002) 


(0 


002) 


F9 


128^ 


24 


-0 


25 





33 





04 





060 





328 





341 





331 


-0 


007 


(-0 


002) 


(0 


003) 


FIO 


128^ 


24 





25 





33 


-0 


04 





060 





342 





328 





329 





007 


(-0 


002) 


(0 


002) 


Fll 


128^ 


24 





50 





33 


-0 


08 





060 





348 





321 





330 





016 


(-0 


001) 


(0 


002) 


F12 


128^ 


24 
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Table B.3. Summary of the runs in the Sets G and H with = -k/2. 
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0.332 


0.342 
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